L’objectif de ce TP est d’illustrer les notions abordées dans le chapitre de régression linéaire.

Les librairies R nécessaires pour ce TP :

library(ellipse)
library(leaps)
library(MASS)
library(corrplot)
library(glmnet)
library(coefplot)
library(ggplot2)  
library(gridExtra)
library(ggfortify)
library(plotly)   
library(reshape2)

1 Introduction

La pollution de l’air constitue actuellement une des préoccupations majeures de santé publique. De nombreuses études épidémiologiques ont permis de mettre en évidence l’influence sur la santé de certains composés chimiques comme le dioxyde souffre (SO2), le dioxyde d’azote (NO2), l’ozone (O3) ou des particules en suspension. Des associations de surveillance de la qualité de l’air (Air Breizh en Bretagne depuis 1994) existent sur tout le territoire métropolitain et mesurent la concentration des polluants. Elles enregistrent également les conditions météorologiques comme la température, la nébulosité, le vent, les chutes de pluie en relation avec les services de Météo France… L’une des missions de ces associations est de construire des modèles de prévision de la concentration en ozone du lendemain à partir des données disponibles du jour : observations et prévisions de Météo France. Plus précisément, il s’agit d’anticiper l’occurrence ou non d’un dépassement légal du pic d’ozone (180 \(\mu\)gr/m3) le lendemain afin d’aider les services préfectoraux à prendre les décisions nécessaires de prévention : confinement des personnes à risque, limitation du trafic routier. Plus modestement, l’objectif de cette étude est de mettre en évidence l’influence de certains paramètres sur la concentration d’ozone (en \(\mu\)gr/m3) et différentes variables observées ou leur prévision. Les 112 données étudiées ont été recueillies à Rennes durant l’été 2001.

Les 13 variables observées sont :

  • maxO3 : Maximum de concentration d’ozone observé sur la journée en \(\mu\)gr/m3
  • T9, T12, T15 : Température observée à 9, 12 et 15h
  • Ne9, Ne12, Ne15 : Nébulosité observée à 9, 12 et 15h
  • Vx9, Vx12, Vx15 : Composante E-O du vent à 9, 12 et 15h
  • maxO3v : Teneur maximum en ozone observée la veille
  • vent : orientation du vent à 12h
  • pluie : occurrence ou non de précipitations

Les données sont disponibles sur la page moodle du cours. Pour les importer, vous pouvez utiliser la commande suivante :

ozone = read.table("Ozone.txt")

Afin de vous familiariser avec les données, faites dans un premier temps une analyse de statistique descriptive. Vous pouvez utiliser les fonctions summary(), boxplot(), pairs(), barplot(), corrplot(), ….

# A COMPLETER - Faire des stat descriptives des données
summary(ozone)
     maxO3              T9             T12             T15       
 Min.   : 42.00   Min.   :11.30   Min.   :14.00   Min.   :14.90  
 1st Qu.: 70.75   1st Qu.:16.20   1st Qu.:18.60   1st Qu.:19.27  
 Median : 81.50   Median :17.80   Median :20.55   Median :22.05  
 Mean   : 90.30   Mean   :18.36   Mean   :21.53   Mean   :22.63  
 3rd Qu.:106.00   3rd Qu.:19.93   3rd Qu.:23.55   3rd Qu.:25.40  
      Ne9             Ne12            Ne15           Vx9         
 Min.   :0.000   Min.   :0.000   Min.   :0.00   Min.   :-7.8785  
 1st Qu.:3.000   1st Qu.:4.000   1st Qu.:3.00   1st Qu.:-3.2765  
 Median :6.000   Median :5.000   Median :5.00   Median :-0.8660  
 Mean   :4.929   Mean   :5.018   Mean   :4.83   Mean   :-1.2143  
 3rd Qu.:7.000   3rd Qu.:7.000   3rd Qu.:7.00   3rd Qu.: 0.6946  
      Vx12             Vx15            maxO3v           vent          
 Min.   :-7.878   Min.   :-9.000   Min.   : 42.00   Length:112        
 1st Qu.:-3.565   1st Qu.:-3.939   1st Qu.: 71.00   Class :character  
 Median :-1.879   Median :-1.550   Median : 82.50   Mode  :character  
 Mean   :-1.611   Mean   :-1.691   Mean   : 90.57                     
 3rd Qu.: 0.000   3rd Qu.: 0.000   3rd Qu.:106.00                     
    pluie          
 Length:112        
 Class :character  
 Mode  :character  
                   
                   
 [ getOption("max.print") est atteint -- ligne 1 omise ]

2 Régression linéaire simple

Dans cette section, nous souhaitons expliquer la concentration d’ozone maximale de la journée (maxO3) en fonction de la température à 12h (T12).

  1. Représentez le nuage de points \((x_i,y_i)\) à l’aide de la fonction plot() (ou geom_point() de ggplot2).
# A completer
# Scatter plot
plot(ozone$T12, ozone$maxO3, 
     xlab = "Temperature at 12:00 PM (T12)",
     ylab = "Maximum Ozone Concentration (maxO3)",
     main = "Scatter Plot of maxO3 vs. T12")

Question : Pensez-vous que l’ajustement d’un modèle de régression linéaire \(y_i = \theta_0+\theta_1 x_i +\varepsilon_i\) est justifié graphiquement?

Réponse : On voit une espèce de tendance linéaire : un modèle de régression linéaire est justifiée.

  1. Effectuez la régression linéaire à l’aide de la fonction lm() et exploitez les résultats.
# Régression linéaire simple
reg.simple = lm(maxO3 ~ T12, data = ozone)

# Afficher un résumé des résultats de la régression
summary(reg.simple)

Call:
lm(formula = maxO3 ~ T12, data = ozone)

Residuals:
    Min      1Q  Median      3Q     Max 
-38.079 -12.735   0.257  11.003  44.671 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -27.4196     9.0335  -3.035    0.003 ** 
T12           5.4687     0.4125  13.258   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 17.57 on 110 degrees of freedom
Multiple R-squared:  0.6151,    Adjusted R-squared:  0.6116 
F-statistic: 175.8 on 1 and 110 DF,  p-value: < 2.2e-16

Interprétation :

Que contiennent les sorties suivantes :

reg.simple$coefficients  
reg.simple$residuals

Réponse : coefficients -> vecteur des valeurs observées de \(\^{\theta}^{obs}\) residuals -> vecteur des résidus

  1. A l’aide des commandes suivantes, tracez l’estimation de la droite de régression sur le nuage de points ainsi qu’un intervalle de confiance à \(95\%\) de celle-ci :
ggplot(ozone, aes(T12, maxO3))+
    geom_point() +
    geom_smooth(method=lm, se=TRUE)+
    xlab("T12")+  ylab("maxO3")

Ce graphique permet visuellement de vérifier l’ajustement des données au modèle de régression linéaire. Que remarquez-vous?

  1. A l’aide de la commande suivante, étudiez les résidus
autoplot(reg.simple,which=c(1,2,4),label.size=2)     

On prendra soin de comprendre les différents graphiques obtenus.

  1. On s’intéresse maintenant à la qualité de prédiction du modèle. On va donc tracer un intervalle de confiance des prédictions à \(95\%\) avec les commandes suivantes. Commentez.
temp_var = predict(reg.simple, 
                  interval="prediction")
new_df = cbind(ozone, temp_var)
ggplot(new_df, aes(T12, maxO3))+
    geom_point() +
    geom_line(aes(y=lwr), color = "red", 
                       linetype = "dashed")+
    geom_line(aes(y=upr), color = "red", 
                       linetype = "dashed")+
    geom_smooth(method=lm, se=TRUE)+
    xlab("T12")+  
    ylab("maxO3")
  1. On va maintenant s’intéresser à la construction d’intervalles de confiance pour \(\theta_0\) et \(\theta_1\) à \(95\%\).
    A l’aide de la fonction confint(), construisez un intervalle de confiance pour chaque paramètre séparément.
# A COMPLETER
# Intervalles de confiance pour theta_0 (intercept) et theta_1 (pente)
IC = confint(reg.simple)
IC

Pour tenir compte de la dépendance entre \(\theta_0\) et \(\theta_1\), on peut aussi construire une région de confiance pour le vecteur \(\theta=(\theta_0,\theta_1)'\) grâce aux commandes suivantes. Comparez les résultats.

df1 = as.data.frame(
             rbind(coefficients(reg.simple),
             ellipse(reg.simple,level=0.95)))
colnames(df1) = c("intp", "slope")
ggplot(data=df1[-1,],aes(x=intp, y=slope))+
  geom_path()+
  annotate("rect",xmin=IC[1,1],xmax=IC[1,2],
  ymin=IC[2,1],ymax=IC[2,2],fill="red",alpha=0.1)+
  geom_point(data=df1[1,], aes(x=intp, y=slope), 
               pch=3)

3 Régression linéaire multiple

Dans cette partie, nous souhaitons analyser la relation entre le maximum journalier de la concentration d’ozone (maxO3) et

  • la température à différentes heures de la journée (T9, T12, T15),
  • la nébulosité à différentes heures de la journée (Ne9, Ne12, Ne15),
  • la projection du vent sur l’axe Est-Ouest à différentes heures de la journée (Vx9, Vx12,Vx15),
  • la concentration maximale d’ozone de la veille (maxO3v)

On va donc utiliser le sous-jeu de données suivant

ozone1 = ozone[,1:11]

et mettre en place un modèle de régression linéaire multiple.

  1. Faites une analyse descriptive de ce jeu de données.
# Afficher un résumé statistique des variables
summary(ozone1)

# Afficher des boîtes à moustaches pour visualiser la distribution des variables numériques
boxplot(ozone1)

# Afficher une matrice de nuages de points pour explorer les relations entre les variables
pairs(ozone1)
  1. Rappelez l’écriture mathématique du modèle de régression linéaire multiple.

\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + ... + \beta_p X_p + \varepsilon \]

  1. Ajustez un modèle de régression linéaire multiple à l’aide de la commande lm() et commentez les résultats (on appelle reg.mul la sortie de R dans la suite).
reg.mul = lm(maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 + maxO3v, data = ozone1)
resume.mul = summary(reg.mul)
summary(reg.mul)

Call:
lm(formula = maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + 
    Vx12 + Vx15 + maxO3v, data = ozone1)

Residuals:
    Min      1Q  Median      3Q     Max 
-53.566  -8.727  -0.403   7.599  39.458 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.24442   13.47190   0.909   0.3656    
T9          -0.01901    1.12515  -0.017   0.9866    
T12          2.22115    1.43294   1.550   0.1243    
T15          0.55853    1.14464   0.488   0.6266    
Ne9         -2.18909    0.93824  -2.333   0.0216 *  
Ne12        -0.42102    1.36766  -0.308   0.7588    
Ne15         0.18373    1.00279   0.183   0.8550    
Vx9          0.94791    0.91228   1.039   0.3013    
Vx12         0.03120    1.05523   0.030   0.9765    
Vx15         0.41859    0.91568   0.457   0.6486    
maxO3v       0.35198    0.06289   5.597 1.88e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 14.36 on 101 degrees of freedom
Multiple R-squared:  0.7638,    Adjusted R-squared:  0.7405 
F-statistic: 32.67 on 10 and 101 DF,  p-value: < 2.2e-16
  1. Etudiez les résidus
autoplot(reg.mul,which=c(1,2,4),label.size=2)     

4 Sélection de variables

Dans le summary(reg.mul), un test est fait sur chaque coefficient. Il revient à tester que la variable n’apporte pas d’information supplémentaire sachant que toutes les autres variables sont dans le modèle. Il est donc préférable d’utiliser des procédures de choix de modèles pour sélectionner les variables significatives. On va ici comparer la sélection de variable obtenue par différents critères: BIC, AIC, \(R^2\) ajusté, Cp de Mallows. Pour cela, vous pouvez utiliser la fonction regsubsets() de la librairie leaps et la fonction stepAIC() de la librairie MASS. Commentez les résultats obtenus avec les différents critères, vous pourrez vous aider des commandes suivantes :

# Utiliser regsubsets pour sélectionner les variables avec différents critères
choix = regsubsets(maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 + maxO3v, data = ozone1)
plot(choix, scale = 'Cp')

plot(choix, scale = 'adjr2')

plot(choix, scale = 'bic')

# Utiliser stepAIC pour sélectionner le modèle basé sur l'AIC
stepAIC(reg.mul)
Start:  AIC=607.26
maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 + 
    maxO3v

         Df Sum of Sq   RSS    AIC
- T9      1       0.1 20827 605.26
- Vx12    1       0.2 20827 605.26
- Ne15    1       6.9 20834 605.30
- Ne12    1      19.5 20847 605.36
- Vx15    1      43.1 20870 605.49
- T15     1      49.1 20876 605.52
- Vx9     1     222.6 21050 606.45
<none>                20827 607.26
- T12     1     495.5 21323 607.89
- Ne9     1    1122.6 21950 611.14
- maxO3v  1    6459.6 27287 635.51

Step:  AIC=605.26
maxO3 ~ T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 + maxO3v

         Df Sum of Sq   RSS    AIC
- Vx12    1       0.1 20827 603.26
- Ne15    1       6.9 20834 603.30
- Ne12    1      22.1 20849 603.38
- Vx15    1      43.5 20871 603.49
- T15     1      49.4 20877 603.52
- Vx9     1     264.0 21091 604.67
<none>                20827 605.26
- T12     1     641.4 21469 606.66
- Ne9     1    1183.6 22011 609.45
- maxO3v  1    7112.2 27940 636.16

Step:  AIC=603.26
maxO3 ~ T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx15 + maxO3v

         Df Sum of Sq   RSS    AIC
- Ne15    1       6.9 20834 601.30
- Ne12    1      23.1 20850 601.38
- T15     1      50.0 20877 601.53
- Vx15    1      79.6 20907 601.69
- Vx9     1     320.1 21148 602.97
<none>                20827 603.26
- T12     1     647.3 21475 604.69
- Ne9     1    1185.1 22013 607.46
- maxO3v  1    7112.6 27940 634.16

Step:  AIC=601.3
maxO3 ~ T12 + T15 + Ne9 + Ne12 + Vx9 + Vx15 + maxO3v

         Df Sum of Sq   RSS    AIC
- Ne12    1      16.2 20851 599.38
- T15     1      45.2 20880 599.54
- Vx15    1      75.1 20910 599.70
- Vx9     1     325.2 21160 601.03
<none>                20834 601.30
- T12     1     971.8 21806 604.40
- Ne9     1    1215.8 22050 605.65
- maxO3v  1    7106.3 27941 632.17

Step:  AIC=599.38
maxO3 ~ T12 + T15 + Ne9 + Vx9 + Vx15 + maxO3v

         Df Sum of Sq   RSS    AIC
- T15     1      45.6 20896 597.63
- Vx15    1      77.6 20928 597.80
- Vx9     1     337.6 21188 599.18
<none>                20851 599.38
- T12     1    1034.0 21885 602.80
- Ne9     1    2155.5 23006 608.40
- maxO3v  1    7130.2 27981 630.33

Step:  AIC=597.63
maxO3 ~ T12 + Ne9 + Vx9 + Vx15 + maxO3v

         Df Sum of Sq   RSS    AIC
- Vx15    1      74.1 20970 596.02
- Vx9     1     366.5 21263 597.58
<none>                20896 597.63
- Ne9     1    2241.1 23137 607.04
- T12     1    6712.6 27609 626.83
- maxO3v  1    7381.5 28278 629.51

Step:  AIC=596.02
maxO3 ~ T12 + Ne9 + Vx9 + maxO3v

         Df Sum of Sq   RSS    AIC
<none>                20970 596.02
- Vx9     1     903.4 21874 598.75
- Ne9     1    2714.8 23685 607.66
- T12     1    6650.4 27621 624.88
- maxO3v  1    7363.5 28334 627.73

Call:
lm(formula = maxO3 ~ T12 + Ne9 + Vx9 + maxO3v, data = ozone1)

Coefficients:
(Intercept)          T12          Ne9          Vx9       maxO3v  
    12.6313       2.7641      -2.5154       1.2929       0.3548  

Avec le critère BIC, nous retenons 4 variables : T12, Ne9, Vx9 et maxO3v.
A l’aide des commandes suivantes, testez le sous-modèle avec les 4 variables retenues par BIC contre le modèle complet. Commentez.

reg.fin=lm(maxO3~T12+Ne9+Vx9+maxO3v,data=ozone1)
summary(reg.fin)
anova(reg.fin,reg.mul)

5 Régressions régularisées

On commence par centrer et réduire les données avant de mettre en place et comparer des méthodes de régression régularisées.

tildeY=scale(ozone1[,1],center=T,scale=T)
tildeX=scale(ozone1[,-1],center=T,scale=T)

5.1 Régression Ridge

  1. A l’aide de la fonction glmnet(), ajustez une régression ridge en faisant varier \(\lambda\) sur une grille. On stockera le résultat dans la variable fitridge. Explorez le contenu de fitridge.
lambda_seq<-10^(seq(-4,4,0.01))

# Ajuster la régression ridge avec glmnet
fitridge <- glmnet(x = tildeX, y = tildeY, alpha = 0, lambda = lambda_seq, family = c("gaussian"), intercept = F)

summary(fitridge)
          Length Class     Mode   
a0         801   -none-    numeric
beta      8010   dgCMatrix S4     
df         801   -none-    numeric
dim          2   -none-    numeric
lambda     801   -none-    numeric
dev.ratio  801   -none-    numeric
nulldev      1   -none-    numeric
npasses      1   -none-    numeric
jerr         1   -none-    numeric
offset       1   -none-    logical
call         7   -none-    call   
nobs         1   -none-    numeric
  1. Tracez les chemins de régularisation de chaque variable et commentez.
df=data.frame(lambda = rep(fitridge$lambda,ncol(tildeX)), theta=as.vector(t(fitridge$beta)),variable=rep(colnames(tildeX),each=length(fitridge$lambda)))
g1 = ggplot(df,aes(x=lambda,y=theta,col=variable))+
  geom_line()+
  theme(legend.position="bottom")+
  scale_x_log10()
ggplotly(g1)
  1. A l’aide de la fonction cv.glmnet() mettez en place une validation croisée pour sélectionner le “meilleur” \(\lambda\) par MSE.
# Réaliser une validation croisée pour sélectionner le meilleur lambda
ridge_cv <- cv.glmnet(x = tildeX, y = tildeY, alpha = 0, lambda = lambda_seq, nfolds = 10, type.measure = c("mse"), intercept = F)

# Trouver la meilleure valeur de lambda en minimisant le MSE
best_lambda <- ridge_cv$lambda.min
df2=data.frame(lambda=ridge_cv$lambda,MSE=ridge_cv$cvm,cvup=ridge_cv$cvup,cvlo=ridge_cv$cvlo)
gmse<-ggplot(df2)+
  geom_line(aes(x=lambda,y=MSE))+
  geom_vline(xintercept = ridge_cv$lambda.min,col="red",linetype="dotted")+
  geom_line(aes(x=lambda,y=cvup),col="blue",linetype="dotted")+
  geom_line(aes(x=lambda,y=cvlo),col="blue",linetype="dotted")+
  #xlim(c(0,ridge_cv$lambda.min+0.5))+
  scale_x_log10()
ggplotly(gmse)

La valeur de \(\lambda\) sélectionnée vaut

g1=g1 + 
  geom_vline(xintercept = best_lambda,linetype="dotted", color = "red")+
  scale_x_log10()

5.2 Régression Lasso

  1. A l’aide de la fonction glmnet(), ajustez une régression Lasso en faisant varier \(\lambda\) sur une grille. On stockera le résultat dans la variable fitlasso. Explorez le contenu de fitlasso.
# Ajuster la régression Lasso avec glmnet
fitlasso <- glmnet(x = tildeX, y = tildeY, alpha = , lambda = lambda_seq, family = c("gaussian"), intercept = F)

# Afficher un résumé des résultats
summary(fitlasso)
          Length Class     Mode   
a0         801   -none-    numeric
beta      8010   dgCMatrix S4     
df         801   -none-    numeric
dim          2   -none-    numeric
lambda     801   -none-    numeric
dev.ratio  801   -none-    numeric
nulldev      1   -none-    numeric
npasses      1   -none-    numeric
jerr         1   -none-    numeric
offset       1   -none-    logical
call         6   -none-    call   
nobs         1   -none-    numeric
  1. Tracez le chemin de régularisation de chacune des variables et commentez
df=data.frame(lambda = rep(fitlasso$lambda,ncol(tildeX)), theta=as.vector(t(fitlasso$beta)),variable=rep(colnames(tildeX),each=length(fitlasso$lambda)))
g3 = ggplot(df,aes(x=lambda,y=theta,col=variable))+
  geom_line()+
  theme(legend.position="bottom")+
  scale_x_log10()
ggplotly(g3)
  1. A l’aide de la fonction cv.glmnet() mettez en place une validation croisée pour sélectionner le “meilleur” \(\lambda\) par MSE. En pratique, il est préconisé d’utilisé lambda.1se (la plus grande valeur de \(\lambda\) telle que l’erreur standard se situe à moins de 1 de celle du minimum).
lasso_cv <- cv.glmnet(...) # A COMPLETER
best_lambda <-lasso_cv$lambda.min
lambda1se <- lasso_cv$lambda.1se

La valeur de \(\lambda\) sélectionnée vaut ….

g3=g3 + 
  geom_vline(xintercept = best_lambda,linetype="dotted", color = "red")+
  geom_vline(xintercept = lambda1se,linetype="dotted", color = "blue")+
  scale_x_log10()
g3
  1. Quelle sélection de variables obtient-on alors ?

Vous pouvez utiliser la fonction extract.coef() de la librairie coefplot

# A COMPLETER

5.3 Régression Elastic Net

  1. Reprenez les questions précédentes pour ajuster une régression Elastic Net
fitEN <- glmnet(...)
df=data.frame(lambda = rep(fitEN$lambda,ncol(tildeX)), theta=as.vector(t(fitEN$beta)),variable=rep(c(colnames(tildeX)),each=length(fitEN$lambda)))
g4 = ggplot(df,aes(x=lambda,y=theta,col=variable))+
  geom_line()+
  theme(legend.position="bottom")+
  scale_x_log10()
EN_cv <- cv.glmnet(...) # A COMPLETER
best_lambda <-EN_cv$lambda.min
g4=g4 + geom_vline(xintercept = best_lambda,linetype="dotted", 
                color = "red")
ggplotly(g4)
  1. Comparez les coefficients obtenus avec la régression linéaire, la régression ridge, la régression Lasso et la régression ElasticNet